1. IntroductionThe CO2 output produced by human activities has a great impact on the Earth’s environment. During the past two centuries, the atmospheric CO2 concentration has been increased by approximately 40% due to fossil fuel combustion, automobile emission and deforestation.[1] About a third of anthropogenic carbon added to the atmosphere was absorbed by the ocean, which greatly affects the carbonate chemistry in the ocean.[2] One of the main reactions for CO2 dissolving in water is: 
. In addition, this reaction also plays a significant role in the pH regulation of blood and transportation of CO2 in living systems.[3] An interesting aspect about the title reaction is that it does not have a transition state in the gas phase,[4–12] however, it does need to overcome a reaction barrier[4,8–13] in an aqueous solution. Namely, the solvent induces a transition state for this reaction in the aqueous solution. Therefore, understanding the mechanism of the water-induced transition state on the atomic level will expand our knowledge of this type of reaction wherein a solvent alters the reaction mechanism by inducing a transition state on the potential energy surface.
Liang and Lipscomb’s ab initio study[7] of the title reaction in the gas phase suggests that the appearance of the barrier in the solution may be associated with displacement of water from hydration stabilized 
. Nemukhin et al. studied this reaction in an aqueous solution using a continuum solvation model and a cluster model consisting of 30 water molecules with an effective fragment potential (EFP) approach,[4] they found that the hydroxide anion losing water molecules from its hydration shell results in the potential barrier, and they estimated the barrier height within the limits 8 kcal/mol–13 kcal/mol with the LMP2 and B3LYP approximations. In addition, various theoretical calculations including variants of continuum models of solvation,[4,6,12] molecular dynamics simulations[8–10] and a RISM-SCF method[11] calculated the free energy activation barrier with its value ranging from 6.1 kcal/mol to 19.2 kcal/mol. However, according to the summary of experimental data,[14] the free energy reaction barrier for this reaction in an aqueous solution is deduced at 12.1 kcal/mol.[4] Therefore, till now, different theoretical studies have produced a big range of uncertainty on the free energy reaction barrier compared to the experimental value.
Due to previous studies mainly using continuum solvation models and a cluster model only consisting of 30 water molecules to describe the aqueous solution,[4,6,12] in this work, we used an explicit, SPC/E water model[15] with 1365 water molecules to treat the aqueous solution to understand the solvent-induced reaction barrier on the atomic level. Furthermore, due to the big range uncertainty, from 6.1 kcal/mol to 19.2 kcal/mol, of the free energy reaction barrier from previous studies,[4–13] here we employed the multi-level, quantum mechanics (QM) theories to treat the solute region to obtain an accurate prediction on the free energy reaction barrier with the CCSD(T)[16,17] level of theory. The multi-level quantum mechanics theories, up to the “gold-standard” CCSD(T) theory, were employed to treat the solute region, where the electrostatic potential (ESP),[18] density function theory (DFT)[19,20] and CCSD(T) levels of theory were used at different stages of the calculation to obtain an accurate potential of mean force (PMF) efficiently. Thus the combined multi-level quantum mechanics theories and molecular mechanics (ML-QM/MM)[21,22] with 1365 SPC/E water molecules was employed to study the 
reaction in an aqueous solution: the accurate PMF along the reaction pathway was mapped at the CSSD(T)/MM level of theory; the detailed, atomic level, solvation-induced reaction mechanism was revealed; and the solvent effect contributions on the PMF and on the solvation-induced transition state were investigated.
2. Methods2.1. ML-QM/MM methodFor the reaction system 
in aqueous, the 
+CO2 is treated as the QM region, and the aqueous is treated as the MM region. Due to the formidable computational cost, it is not realistic to finish the whole dynamical calculation process under the CCSD(T) level of theory. Hence a multi-level QM theory, including ESP,[21] DFT, and CCSD(T), was used to treat the QM region during the different stages of our QM/MM calculation to eventually obtain the more accurate PMF with the “gold-standard” CCSD(T) level of theory. Consequently, the ML-QM/MM approach including CCSD(T)/MM, DFT/MM, and ESP/MM is utilized to treat the whole reaction system. This ML-QM/MM approach has been proven to achieve more accurate free energy barrier heights for reactions in the solution phase[23,24] than the usual DFT treatment on the QM part.
The potential of mean force along the reaction path for two adjacent configurations A and B under the CCSD(T)/MM level of theory can be expressed as
The first term 
is the MM contribution to the PMF at the ESP/MM level of theory at fixed A and B solute configurations. The second term and third term in brackets denote free energy differences for shifting from the ESP to DFT, and from DFT to CCSD(T) levels of theory.
2.2. Numerical detailsThe reactant (
) was embedded into a 34.6-Å cubic box of 1365 water molecules described by an explicit, SPC/E water model. The QM region was computed with DFT and CCSD(T) levels of theory at different stages of calculation. The aug-cc-pvDZ basis set was used for both DFT and CCSD(T) calculations, and the M06-2X exchange correlation function was used for the DFT theory. A cutoff radius of 15 Å was used to separate the bonded, electrostatic and Van der Waals interactions between the solute and solvent. The van der Waals parameters for the QM region were taken from a standard Amber force field.[25]
The structure of the initial reactant complex was set up based on the structure optimized at the MP2/6-311++G** level by Peng et al.[9] and was embedded into the water box described above. The whole reaction system was optimized with a multi-region optimization approach. After this, the water molecules (MM region) were equilibrated using the molecular dynamics simulation for 40 ps with a time step of 0.001 ps at 298 K, while the reactants (QM region) were fixed and represented by the ESP charges obtained from the previous step. The initial reactant complex structure in water was obtained after the whole system was optimized again. Following the 
bond forming reaction mechanism, we obtained the initial product complex, which was also equilibrated and optimized to reach its equilibrium.
Next, the initial reaction pathway was constructed using the climbing image NEB method with the initial reactant and product complexes obtained above. The highest energy point along the initial reaction pathway was isolated for saddle point search. The transition state was found and confirmed by a numerical frequency calculation with only one imaginary frequency at 
. The vibrational normal mode of the imaginary frequency was displaced toward the reactant and product sides respectively, and the multi-region optimization was performed on the displacements again to obtain the final reactant and product complexes.
Finally, the reaction pathway was generated with the final reactant and product complexes with 10-bead NEB calculations. The solvent along the pathway for each bead was optimized and the molecular dynamics simulation was carried out on each NEB bead for 40 ps, then optimized again. The above steps were repeated until we obtained the converged reaction path. Then the PMFs with DFT/MM and CCSD(T)/MM levels of theory were calculated according to Eq. (1).
2.3. Gas phase reaction pathIn order to verify our results on the CCSD(T)/MM level of theory are reliable and accurate, we also used the gas-phase reaction path and Born solvation model to calculate a reaction profile in solution to compare with our CCSD(T)/MM reaction path. The gas-phase reactant state was constructed by adopting the experimental values[26,27] for its 
and CO2 structures and they were set apart at 10.0 Å; the structure of the last bead on the NEB reaction path in the aqueous solution was optimized in the gas phase at the DFT/M06-2X/aug-cc-pVDZ level of theory as the gas-phase product complex. Then the gas-phase reaction path was generated with the above reactant state and product complex using the NEB method. Ten NEB beads in the gas phase, which correspond to the 10 NEB beads in the aqueous solution, were obtained first using the above DFT level of theory; then we shifted to the CCSD(T) domain to calculate the reaction path energy using the CCSD(T)/aug-cc-pVDZ level of theory.
3. Results and discussion3.1. Potential of mean force and solvent energy contributionThe PMFs of the reaction 
in water along the nudged elastic band (NEB[28]) reaction path with both the DFT/MM and CCSD(T)/MM levels of theory, as well as the solvent contribution and the internal energies, are plotted in Fig. 1. Note, only relative free energy along the reaction path is calculated by taking the reactant free energy as a zero reference energy point. First, there is already no energy barrier for the internal energies(black and red lines) without the solvent energy contribution; however, there is a solvent-energy-induced energy barrier appearing at the No. 7 bead in the aqueous solution for either level of theory, DFT/MM or CCSD(T)/MM. The numerical frequency calculation for this point in the aqueous solution showed a single imaginary frequency of 
, thus the No. 7 bead is indeed the induced transition state for this reaction in the aqueous solution. It is not a coincidence that there is a maximum of the solvent energy contribution, 42.8 kcal/mol, exactly at the seventh point along the reaction path, causing the No. 7 bead to have the highest point on the PMFs. For example, for the CCSD(T)/MM level of theory, this maximum contribution, 42.8 kcal/mol, outweighs the energy decrease of the internal energy, −31.0 kcal/mol, at the No. 7 bead, thus leading to the highest energy point on the PMF. Therefore, in terms of energy, it is solely the solvent energy contribution that caused the transition state to be the highest point on the PMF.
Second, the free energy barrier in the solution for this reaction at the DFT/MM level of theory is 8.8 kcal/mol, it reaches 11.8 kcal/mol when shifted from the DFT/MM to the CCSD(T)/MM level of theory. The free energy barrier height at the CCSD(T)/MM level of theory is in very good agreement with the experimental value at 12.1 kcal/mol.[4,14] The reaction free energy of this reaction in water is −25.1 kcal/mol at the DFT/MM level of theory and shifted to −17.7 kcal/mol at the CCSD(T)/MM level of theory, the latter is very close to the one around −15.2 kcal/mol estimated by Nemukhin’ group[4] using a continuum model with DFT/B3LYP/6-311++G(d, p) theory.
3.2. Comparison of the reaction profilesIn order to be consistent with the PMF in aqueous solution, the NEB reaction path in the gas phase was also constructed using a multi-level quantum mechanics method with DFT/M06-2X/aug-cc-pVDZ and CCSD(T)/aug-cc-pVDZ levels of theory. Figure 2 shows the schematic plot of the comparison of the NEB reaction paths between the gas phase(with CCSD(T) theory) and aqueous solution(with CCSD(T)/MM theory), as well as the reaction stationary points in aqueous solution obtained using free energies of solvation from the gas phase. The calculated gas-phase reaction is an exoergic reaction with a reaction energy at −46.3 kcal/mol, and it indeed does not have a transition state as the reaction path goes downhill from bead 1 to bead 10. The energy difference between the No. 7 bead and the reaction state in the gas phase is −23.5 kcal/mol under the CCSD(T) level of theory. This Figure also indicates, using the gas-phase reaction path as the starting point, that it is the solvation energy difference among the stationary points that results in the energy barrier for this reaction in solution.
3.3. From gas phase to solution phaseFor the gas phase reactant state and product state, 
, CO2, and 
, 
is the most compact molecule with charge among the three, thus it is the easiest to be solvated in an aqueous solution. The estimated data from experiments[9] shows that the 
has the biggest free energy of solvation at −113 kcal/mol, then 
at −81.5 kcal/mol, and the CO2 has the smallest at 0.11 kcal/mol since it is neutral. For the 
complex at the No. 7 bead in the gas phase, which corresponds to the transition state in the aqueous solution, we used the Born solvation model to estimate its free energy of solvation,
where
ε0,
,
q, and
R denote the relative dielectric constant in a vacuum and in water, charge of solute, and the solvation radius. Here, the
Rvalue, 2.128 Å, is the farthest distance from atoms in
complex to its center of mass. With values of
ε0,
,
q are 1.0, 78.5, and −1.0 respectively, the free energy of solvation using this model was calculated at −76.5 kcal/mol for the No. 7 bead. Its absolute value is the smallest among the stationary points on the reaction path and is expected since the No. 7 bead is the least compact molecule among the reactants,
and the product complex
.
The free energies of solvation of 
and CO2 are −113 kcal/mol and 0.11 kcal/mol[9] respectively, which combines to give the free energy of solvation of the reaction at −112.89 kcal/mol (as labeled in Fig. 2); the free energy of solvation of the transition state calculated using Eq. (2) mentioned above is −76.5 kcal/mol; the gas-phase energy of the No. 7 bead is −23.5 kcal/mol; therefore, combining these three numbers above, we deduced the highest energy point at 12.9 kcal/mol 
on the reaction path in aqueous solution (as shown in Fig. 2), which is the transition state of this reaction in the aqueous solution. This estimated transition state height, 12.9 kcal/mol, agrees very well with our ML-QM/MM result at 11.8 kcal/mol at the CCSD(T)/MM level of theory, and the estimated free energy of a reaction using the same process at −14.9 kcal/mol also agrees with our ML-QM/MM result at −17.7 kcal/mol at the CCSD(T)/MM level of theory. More importantly, our calculated energy barrier height at 11.8 kcal/mol at the CCSD(T)/MM level of theory has a very good agreement with the experimental value at 12.1 kcal/mol.
3.4. Evolution of the water induced transition stateThe water-induced transition state appearing in the aqueous solution is not only the change of the energy but also the change of structures due to the solvent-solute interactions. Figure 3 shows the geometry evolution of the 10 beads along the NEB reaction pathway in the aqueous solution. The first bead, the seventh, and the last bead along the NEB reaction pathway correspond to the reaction complex, transition state, and product complex respectively. For the reactant complex, there are five hydrogen bonds formed with the 
with an average 
distance at 1.654 Å; however, there are no hydrogen bonds formed with CO2. Usually the solvation pattern of the 
is with four or five hydrogen bonds[29,30] and CO2 with four hydrogen bonds.[31] Hence, at the reactant complex, the 
forms its solvation shell while CO2 does not. This also means that 
is very well solvated in water, and CO2 is not, which has been proven with the free energies of solvation of 
at −113 kcal/mol and CO2 at 0.11 kcal/mol9. As the reaction proceeds, the 
gradually loses it hydrogen bonds while the CO2 gains hydrogen bonds. This is understandable since CO2 gains a partial negative charge as 
loses its partial negative charge along the reaction path. This can be seen from Fig. 4, as the negative charge of 
decreases, the negative charge of CO2 increases concertedly from bead 1 to bead 7 as 
decreases from −1.0 to −0.75, and CO2 increases from 0.0 to −0.25.
At the transition state No. 7, the strong solvated 
only has four strong hydrogen bonds with an average 
distance at 1.707 Å, while CO2 gains three strong hydrogen bonds with an average distance at 1.791 Å. Thus, from reactant state (Fig. 3(a)) to transition state (Fig. 3(g)), the 
still has its solvation shell while CO2 gains hydrogen bonds however without forming a complete solvation shell. Then the largest charge change occurs from bead 7 to bead 8, the 
negative charge drops significantly from −0.75 to −0.37 and at the same time CO2 has gained the most charge at this step from −0.25 to −0.63. As a result, the 
loses its solvation shell with only one hydrogen bond left, while CO2 forms its solvation shell with four hydrogen bonds with an average distance at 1.694 Å. Hence No. 7 bead serves as the turning point of one solvation shell (OH−) breaking and the other solvation shell (CO2) formation: at No. 8 bead, 
has lost three hydrogen bonds to destroy its solvation shell while CO2 has added another hydrogen bond to form its solvation shell; the loss of the hydrogen bonds of OH− cannot be compensated by the formation of the new hydrogen bonds of CO2. Thus the solvent energy of the No. 7 bead is the largest on the solvent energy contribution curve as shown in Fig. 1, which induced an energy barrier in solution due to the solvation shell changes. Till the end of the reaction, 
holds its one water molecule while CO2 has its new formed solvation shell with four water molecules. Note another feature of the transition state is the 
position, there is an inversion of the OH− position before and after the transition state: before the transition state, 
titles towards one of the O atom in CO2; at the transition state, 
is almost perpendicular to the O (as in 
distance with the 
angle at 92.9°; after the transition state, the 
direction is inverted to title away from the strong-bended CO2 with an 
angle at 129.8°.
4. ConclusionThe solvent energy contribution has a maximum on its reaction path, and this maximum induces the highest point, the energy barrier, on the reaction path in an aqueous solution. From the structure point of view, one has to take into account of the solvation shell changes of the two reactants: the broken of the 
solvation shell and the formation of the CO2 solvation shell during the reaction process induces the transition state in the aqueous solution. The induced barrier height from our calculation at the CCSD(T)/MM level of the theory is 11.8 kcal/mol, which agrees very well with the experimental result, 12.1 kcal/mol. For the title reaction in solution, the reaction mechanism is altered by the presence of the aqueous solution due to solvation-shell changes.
